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ABSTRACT. Solutions of a diophantine equation f(a,b) = g(c,d), with a,b,c,d in some 
finite range, can be efficiently enumerated by sorting the values of / and g in ascending 
order and searching for collisions. This article considers functions /: MxH^Z that 
are bimonotone in the sense that f(a,b) < f(a',b') whenever a < a' and b < b'. A two- 
variable polynomial with non-negative coefficients is a typical example. The problem is 
to efficiently enumerate all pairs (a,b) such that the values f(a,b) appear in increasing 
order. We present an algorithm that is memory-efficient and highly parallelizable. In 
order to enumerate the first n values of /, the algorithm only builds up a priority queue of 
length at most \f%n + 1 . In terms of bit-complexity this ensures that the algorithm takes 
time 0(nlog 2 n) and requires memory 0(\/n\ogn), which considerably improves on the 
memory bound @(«log«) provided by a naive approach, and extends the semimonotone 
enumeration algorithm previously considered by R.L. Ekl and D.J. Bernstein. 



1. Introduction and statement of results 

1.1. Motivation. Given polynomial functions f,g: N x N — > Z, how can we efficiently 
enumerate solutions of the equation f(a,b) = g(c,d)l One standard way to do this is to 
sort the sets F = { (f(a,b),a,b) | 1 < a,b < N} andG = { (g(c,d),c,d) | 1 <c,d <jV}into 
ascending order with respect to the first coordinate and to look for collisions. As stated, this 
requires to store all elements before sorting, which consumes memory 0(nlog«), where 
n =N 2 is the number of values to enumerate, and time between £2(nlogn) and 0(nlog 2 n). 

The present article develops a less memory consuming algorithm under the hypothesis 
that / and g are bimonotone, that is, monotone in each variable. This is sufficiently of- 
ten the case to be interesting, for example when / and g are given by polynomials with 
non-negative coefficients. Given a bimonotone function /, Algorithm|U discussed below, 
produces a stream Jti,JC2,Jt3j ■ ■ ■ enumerating all parameters Xj — (a,,/?,) in the domain of 

/ such that f(x\) ^ f{xi) ^ f( x 3) ^ Having at hand such sorted enumerations for 

/ and g, one can easily enumerate solutions of the equation f(x) = g(y): start with ; = 1 
and 7 = 1; whenever /(x,) < giyj), increment i; whenever /(x,) > g(yj), increment /'. If 
eventually /(x,) = giyi), then output the solution (x^yA and continue searching. 

1.2. Main result. The idea of sorted enumeration has been applied by D.J. Bernstein JT] 
with great success to equations of the special form p(a)+q(b) — r(c) +s(d). We generalize 
his approach to arbitrary bimonotone functions. The main result can be stated as follows: 



Date: December 19, 2009. 

2000 Mathematics Subject Classification. 68P10; 11Y50, 68W10, 11Y16, 11D45. 

Key words and phrases, sorting and searching, diophantine equation, bimonotone function, sorted enumera- 
tion, semimonotone enumeration, bimonotone enumeration. 

1 



2 



MICHAEL EISERMANN 



Theorem 1. Suppose that f: NxM->Zis bimonotone and proper in the sense that for 
every z G Z only finitely many pairs (a,b) satisfy f(a,b) < z. Then Algorithm^ stated 
below produces a stream enumerating all pairs (a,b) € N x N such that the values f{a,b) 
appear in increasing order. While enumerating the first n values, the algorithm only builds 
up a priority queue of length m < \/2n + 1. If f is a polynomial, this ensures that the 
algorithm takes time 0(n\og 2 n) and requires memory 0(^/n\ogn). 

The precise bound m < \fln + 1 is free of hidden constants and thus uniformly valid 
for all bimonotone functions /. The less explicit bounds of time 0(n\og 2 n) and memory 
O(yfnlogn) concern the bit-complexity and the hidden constants necessarily depend on /. 
We shall assume throughout that / behaves polynomially, see £12.31 

To place this result into perspective, notice that the time requirement O(n\og 2 n) is 
nearly optimal: enumerating n elements obviously needs n iterations, and one logn factor 
is due to their increasing size. On the other hand, the standard approach would require 
memory ©(nlogn) to store all values before outputting them. Here the stream approach 
can achieve considerable savings and reduce memory to O(^jn\ogn). 

Example 2. Consider f(a,b) — p(a) +q(b) where p and q are non-decreasing polyno- 
mial functions of degree a = degp and j3 = degq, respectively. Assuming 1 < a < j3, 
Algorithm|4]builds up a priority queue of length m 6 &(n e ) with e = 6 [0, j\. 

This illustrates that in the uniform bound m € 0(n l l 2 ), stated in the theorem for all 
bimonotone functions, the exponent i cannot be improved. Notwithstanding, the algorithm 
performs better on certain subclasses of bimonotone functions, where £ <\- 

Remark 3. The predecessor of our algorithm is semimonotone enumeration, recalled in 
S}3] It was devised in [2 1| for polynomials of the form /(a, Z?) = p(a) +q(b), where itpro- 
vides the desired memory bound 0{^/n\ogn). In the more general setting of bimonotone 
functions, however, we show that it only guarantees the memory bound 0(n\ogn) and the 
exponent 1 can in general not be improved. See <2]for a detailed discussion. 

As an additional benefit, our algorithm turns out to be highly parallelizable: 

Remark 4. Algorithm [4] can be adapted to enumerate only those pairs (a,b) G N x N for 
which the values f(a,b) lie in a given interval [zijZa]. Time and memory requirements 
are essentially the same as before; only initialization induces some additional cost and can 
usually be neglected. This means that searching solutions f(a,b) = g(c,d) can be split up 
into disjoint intervals and thus parallelized on independent machines (see fQ. 

1.3. How this article is organized. Section |2] introduces the necessary notation and re- 
calls the generic algorithm of sorted enumeration for an arbitrary map /: X — >Z, where 
X is a finite set. Section [3] discusses a refined algorithm, essentially due to R.L. Ekl J2] 
and D.J. Bernstein |T), under the hypothesis that /: A x B — > Z is semimonotone, that is, 
monotone in the first variable. Section |4] develops a sorted enumeration algorithm for 
bimonotone functions, and Section [5] analyses the asymptotic complexity. Section ^high- 
lights the intrinsically parallel structure of such a search problem. Section [7] generalizes 
our algorithms to functions / : X — > Z restricted to suitable domains X C A x B that are of 
of practical interest. Finally, Section [8] briefly indicates applications to diophantine enu- 
meration problems, such as the taxicab problem. 

1 .4. Acknowledgements. I thank the anonymous referee for his thorough critique, harsh 
but fair, which substantially contributed to improve the exposition. 
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2. Sorted enumeration for arbitrary functions 

Before discussing more sophisticated versions, let us first describe the general problem 
of sorted enumeration and recall its generic solution. 

2.1. The generic problem. Throughout this article we consider an ordered set (Z, ^). By 

order we always mean a reflexive, transitive relation that is complete and antisymmetric, 
i.e. each pair z ^ z! in Z satisfies either z ^ z! or z! ^ Z- Without completeness we may have 
neither z ^ z! nor z! ^ z, in which case we speak of a partial order. Without antisymmetry 
we may have both z ^ z! and z! ^ z, in which case we speak of a preorder. 

We assume that X is a finite or countably infinite set. An enumeration of X is a stream 
x\,X2,x$,... in which each element of X occurs exactly once. Such an enumeration is 

monotone or sorted with respect to /: X — *Z if it satisfies f(x\) ^ fixi) ^ f{x3) ^ 

Whenever the function / is understood from the context, we will simply speak of a sorted 
enumeration ofX. 

Remark 5. The map /: X — > Z can be used to pull back the order ^ from Z to the initially 
unordered set X. More explicitly, we define x =^ x' if and only if f(x) ^ /(*')■ A sorted 
enumeration of X is thus a stream in which all elements of X appear in increasing order 
with respect to the preorder =5;. 

2.2. The generic algorithm. In the general setting, where X is finite and / has no further 
structure, there is essentially only one way to produce a sorted enumeration: 



Algorithm 1 Sorted enumeration for an arbitrary function 
Requires: a function / : X — > Z from a finite set X to an ordered set Z 
Output: an enumeration of X, monotone with respect to / 

1: Generate a list L of all pairs (f(x),x) with x e X. 

2: Sort the list L according to the first coordinate f(x). 

3: Output the arguments x as sorted in the list L. 



Algorithm Q] is obviously correct. Given a set X of size n, generating and reading the 
list L takes n iterations, while sorting requires 0(n\ogn) operations. Not much optimiza- 
tion can be expected concerning these time requirements, since enumeration (sorted or not) 
takes at least n iterations. Memory requirements, however, may be far from optimal, and 
the more specialized algorithms discussed below will mainly be concerned with minimiz- 
ing the use of temporary memory. 

2.3. Time and memory requirements. Throughout this article we use standard asymp- 
totic notation, as in J3] §9]. It is customary to consider the cost for storing and handling 
elements x and f(x) to be constant. This is no longer realistic when the size n = \X\ grows 
without bound. As a typical example, consider a polynomial function / : N — > Z restricted 
to X = {1, . . . ,«}. If each element x G X is stored in binary form, the maximal memory 
required is ©(logn). Likewise, the maximal time to calculate, copy, and compare values 
f(x) is ©(logn), neglecting factors of order log logn or less. Most elements require nearly 
maximum cost, so we shall only consider the worst case. 

In general, we say that / behaves polynomially if the bit-complexity per element is 
0(n), as above. In this case we arrive at the following more realistic account: 

Proposition 6. In order to enumerate a set X of size n, the generic Algorithm\l\builds up 
a list of size m = n, and thus requires time 0(n\og 2 n) and memory of size ©(nlogn). □ 
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3. Sorted enumeration for semimonotone functions 

In this section we consider a semimonotone function /: A x B — > Z. By this we mean 
that (A, ^) is an ordered set and a ^ a' implies f(a,b) ^ f(a',b) for all b £ B. This is the 
same as saying that / is monotone with respect to the partial order (a,b) < (a' ,b') defined 
by the condition a ^ a' and b = b'. 

3.1. The idea of semimonotone enumeration. We will first assume that A and B are 
finite sets. This entails that (A,^) is isotonic to an interval {1,...,/} of integers. The 
minimal and maximal element of A is denoted by and a max , respectively, and the 
successor function is denoted by a aa. Of course a max cannot have a successor in A, so 
by convention we set (7a max = +°°. 

We equip X=Ax5 with the partial order < as defined above. Given a subset X\ C X, 
we denote by M, = MinpQ) the set of its minimal elements. Conversely, M, defines its 
upper set Mf — {x £ X \m <x for some m £ M; }. Figure Q] shows a subset Xj (indicated 
by crosses) together with its set of minima M, (circled crosses). In this example Xj is 
saturated in the sense that Xj — M+. 
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Figure 1 . A subset of A x B and its minima with respect to < 

Since / is monotone with respect to <, the minimum of f(Xi) is attained on M,-. It thus 
suffices to find x\ £ M,- realizing f{x\) — min/(M ( ). We can then output x, and continue 
with the set X{ + \ = X{ \ {x,}. Notice that Xj + \ is again saturated and M,-+i can be easily 
constructed from M,-. This is the idea of Algorithm|2]below. 

3.2. Suitable data structures. The following algorithm has been independently devel- 
oped by R.L. Ekl [2| and D.J.Bernstein [1|, and formalizes the above approach: instead 
of handling the entire set X{, it operates on two smaller sets, M — Min(X,-) and F = f(M). 
In order to efficiently find x, £ M realizing f(xi) = min/(M), we store the set of images 
f(M) in a priority queue F. Recall that a priority queue F for elements of (Z, ^) provides 
the following elementary operations: 

• Inserting an element z £ Z into F ("push"). 

• Reading and removing a minimal element of F ("pop"). 

Priority queues are typically implemented using a heap or a binary tree; in either case 
the elementary operations need 0(logm) steps, where m is the number of elements in the 
priority queue. For a general presentation see Knuth [7 , §5.2.3]. 
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3.3. The semimonotone enumeration algorithm. Instead of /: A x B — > Z it is more 
convenient to work with the map /*: AxB^ZxAxB defined by f* (a, b) = (f(a, b) , a, b). 
In our formulation of Algorithm[2]we thus use a priority queue F for elements inZxAxB, 
sorted by the first coordinate. 



Algorithm 2 Sorted enumeration for a semimonotone function 

Requires: a semimonotone function / : Ax B — > Z 

Output: an enumeration of A x B, monotone with respect to / 

1: Start with an empty priority queue F, then insert /* (a m i n ,b) for all b £ B. 

2: while F is non-empty do 

3: Remove a minimal element f*(a,b) from F and output (a,b). 
4: if era < +oo then insert /* (era, b) into F end if 
5: end while 



Remark 7. All algorithms presented here can be regarded as templates, to be instantiated 
for the given map / : A x B — > Z. Alternatively, one could consider them as taking the sets 
A and B and the map / as input data. In this case, of course, we do not pass the entire sets 
A and B as parameters, nor the map /, say as some subset oiAxBxZ: for finite sets this 
would be as inefficient as AlgorithmQ] for infinite sets it is simply impossible. 

Instead, it suffices to call some function that calculates f(a,b) for any given pair of 
parameters a 6 A and b 6 B. To represent the sets A and B, all we need is the usual iterator 
concept, providing a pointer to the first (and possibly the last) element of the set and a 
method for incrementing, denoted by a i— ► aa above. (Algorithms[5]and[7]also decrement, 
denoted by Z? i — > nb.) 

One can easily add suitable specifications when passing to concrete implementations. 
For the present general exposition, however, we shall maintain the slightly coarser descrip- 
tion, trying to strike a balance between general concepts and implementation details. 

Algorithm[2]is obviously correct. The point is, as motivated above, that it usually uses 
less memory than the generic AlgorithmQ] 

Proposition 8. In order to enumerate a set X — A x B of size n, Algorithm\2\only builds 
up a priority queue of size m = \B\. Let f: M x N ^ Z a semimonotone map that 
behaves polynomially. Applied to subsets A = {I, . . . ,1} and B — {I,..., m} with m > 2, 
the algorithm thus takes time 0{n\ogn\ogm) and requires memory 0(mlog«). 

Proof. The algorithm needs memory to hold m elements f*(a,b) in the priority queue F. 
Since most elements need memory of size ©(logn), we arrive at a total memory cost of 
©(mlogn). During each one of the n iterations, the most time consuming operation is 
updating the priority queue F which requires time O (log n log m). Here m is the size of the 
queue and logn is the typical size of its elements. □ 

Remark 9. Notice that in the degenerate case \B\ = 1, Algorithm [2] simply enumerates 
A in increasing order, which takes time <9(nlogn) and memory ©(logn). In the opposite 
extreme \A\ = 1, it sorts B with respect to / via heap-sort. We thus fall back on the generic 
AlgorithmQ] which takes time (9(mlog 2 m) and space ©(mlogm). 

3.4. Enumerating infinite sets. Sorted enumeration can be generalized from finite to in- 
finite sets. First of all, in order to be amenable to enumeration, A must be either finite 
or isotonic to the natural numbers. Moreover, we have to require that /: A x B — * Z be 
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a proper map in the sense that for every z G Im(/) only finitely many pairs (a,b) satisfy 
f( a ib) < Z- (This condition actually implies thatA is finite or isotonic to N.) Of course, we 
also have to assume that comparisons and all other operations are computable; as before 
their cost will be assumed to be of order (9(logn). 

Proposition 10. Suppose that A is an infinite ordered set and B is a finite set of size m. 
Let f : A x B — > Z be a proper semimonotone map. Then Algorithm\2\produces a stream 
enumerating all pairs (a,b) G A x B such that the values f{a,b) appear in increasing order. 
Producing the first n values takes time 0(n\ogn\ogm) and requires memory ©(wlogn). 

Proof. For every z G Im(/) the set { (a,b) G A x B | f(a,b) < z} is finite, and thus con- 
tained in some finite product [a m i n ,ai] x B. Hence Algorithm |4] correctly enumerates all 
parameters (a,b) with f(a,b) < z as in the finite case. Since this is true for all z, the 
enumeration exhausts Ax B. Bit-complexity behaves as in Proposition^ □ 

We wish to adapt semimonotone enumeration to the case where both A and B are infinite. 
Algorithm|2]is certainly not suited for this task, because the initialization will get stuck in 
an infinite loop. As a necessary restriction we require that / : A x B — > Z be proper, and as 
before we assume that / is monotone with respect to A. For every z G Im(/), we can thus 
enumerate the finite set {/ ^ z} := {(a,b) G A x B \ f(a,b) ^ z} by applying Algorithm^ 
to the relevant finite set B(z) = pr 2 {/ ^ z} = {b G B \ f(a m { n ,b) < z}. 

In order to formulate an explicit algorithm, we assume that the set B is ordered and that 
f2 : B — > Z, /2(b) = f{a m [ n ,b) is non-decreasing. This is strictly weaker than demanding / 
to be bimonotone, because we require monotonicity in b only on the axis {a m in} x B. This 
technical condition ensures that we can easily construct the relevant finite set B(z). In fact, 
the monotonicity of fi is not at all restrictive, because we can choose the order on B, for 
example by pulling back the order on Z via fi to a preorder on B, and then refining to an 
order by arbitrating collisions. In other words, the order on B is just a convenient way to 
encode some preparatory analysis of the proper map fi : B — > Z. 

This idea is formalized in Algorithm [3] which is a slight variation of Algorithm |2] 
The only difference is that it automatically adapts the relevant interval B(z) = [i> m in,^max[ 
according to the level z attained during the enumeration. 



Algorithm 3 Sorted enumeration for a semimonotone function 
Requires: a proper semimonotone function / : Ax B — > Z 
Output: an enumeration of A x B, monotone with respect to /. 

1: Initialize the priority queue F <— (f* (a m j n , b m { n ) ) and set b max <— b t 

2: while F is non-empty do 

3: Remove a minimal element f*(a,b) from F and output (a,b). 
4: if oa < +00 then insert /* (a a, b) into F end if 
5: if b — frmax then 
6: Set 

7: if Vax < +°° then insert /*(<Zmiii,&max) into F end if 
8: end if 
9: end while 



Here we have formulated Algorithm [3] so that it applies to finite and infinite sets alike. 
If A or B is infinite, then aa < +°° or b max < +°°, respectively, is always true and the 
corresponding test can be omitted. 
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Proposition 11. Suppose that both A and B are infinite ordered sets and that f: Ax B —>Z 
is a proper semimonotone function. We also assume that b i— > /(a m in, b) is non-decreasing. 
Then Algorithm\2\provides a sorted enumeration of A x B. While enumerating the first n 
values, it builds up a priority queue of length m < n+ 1. This ensures that the algorithm 
takes time 0(n\og 2 n) and memory 0(n\ogn). □ 

Semimonotone functions are tailor-made for applications where we have monotonicity 
in a but not necessarily in b. They are halfway towards bimonotone functions, which are 
more restrictive but support much better algorithms. These will be discussed next. 

4. Sorted enumeration for bimonotone functions 

In this section we finally turn to bimonotone functions / : A x B — > Z. By this we 
mean that both (A,^) and (B,^) are ordered sets, and that a < a' and b < b' implies 
f(a,b) < f(a ! ,b'). This is the same as saying that / is monotone with respect to the partial 
order (a,b) <^ (a',b') defined by a ^ a' and b ^ b' . 

4. 1 . The idea of bimonotone enumeration. We will first assume that both sets A and B 
are finite. Given a subsetX, C X we denote by M, = Min(X,) the set of its minimal elements 
with respect to <<. Conversely, M, defines its upper set M? = { x £ X \m<<x for some m S 
Mi}. Figure Q] shows a subset Xj (indicated by crosses) together with its set of minima M, 
(circled crosses). In this example X; is saturated in the sense that Xj = Mf. 
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Figure 2. A subset of A x B and its minima with respect to « 



Since / is monotone with respect to <C, the minimum of f(Xj) is attained on M,. It thus 
suffices to find Xi € M; realizing /(x,) = min/(M,). We can then output x, and continue 
with the set Xj + \ = Xi \ {x,}, which is again saturated. Moreover, it is possible to construct 
M; + i directly from M,, without having to construct Xi or Xj+i. (See Algorithm [4] below.) 
Thus, instead of searching the entire set X{, we only need to keep track of M,, the set of 
minimal elements. 

4.2. Suitable data structures. According to the previous remark, the bimonotone enu- 
meration algorithm will operate on two sets: M = Min(X,) and F = f(M). The set M can 
profitably be implemented as a list ( {a\,b\), (a,2,b2),... , (a m ,b m ) ) with at G A and bi £ B. 
During the algorithm, the list M will always be ordered in the sense that a \ < ai < ■ ■ ■ < a m 
and b\ > bi > ■ ■ ■ > b m , as already indicated in Figured We call {a^^bk) the predecessor of 
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{ak+\,bk+\) inM, and conversely (ajt+i,£>jt+i) the successor of (a^b^) inM. By conven- 
tion the predecessor of (ai,b\) is (— °o,+°o), and the successor of (a m ,b m ) is (+°°,— °°). 
Given an element (a,b) in the list M, the required operations are: 

• Finding the successor or predecessor of (a, b) in M. 

• Inserting an element into M right after (a,b). 

• Removing (a,b) from M. 

The cost of these operations can be assumed to be O(\ogn), which is the typical cost 
for storing and handling one of the elements of the set X = A x B of size n. In particular, 
the cost is independent of the size m = \M\. For details on bidirectional lists see Knuth [6„ 
§2.2.5], or any other textbook on algorithms and data structures. 

As before, the set F = f(M) will be implemented as a priority queue containing the 
values f(a,b) for all (a, b) in M. It is recommendable to store f(a,b) together with a 
pointer to the element (a,b) in the list M. This allows us to extract (a,b), and, moreover, 
we can directly address (a,b) in M without searching the list. For notational convenience 
we will not explicitly mention this pointer in the sequel. 

4.3. The bimonotone enumeration algorithm. Having suitable data structures at our dis- 
posal, it is an easy matter to formalize bimonotone enumeration (Algorithm|4]i. 



Algorithm 4 Sorted enumeration for a bimonotone function 

Requires: a bimonotone function / : Ax B — > Z 

Output: an enumeration of A x B, monotone with respect to / 

1: Initialize M «- ( (a m i n ,Vin) ) andF <- ( f{a min ,b min ) >. 

2: while F is non-empty do 

3: Remove a minimal element f(a,b) from F and output (a, b). 
4: Let (a* , b* ) be the successor of (a , b) in the list M. 
5: if aa < a* then 

6: Insert (oa,b) into M right after (a,b) and insert f(oa,b) into F. 
7: end if 

8: Let (a*,b*) be the predecessor of (a,b) in the list M. 
9: if ab < b* then 

10: Insert (a, cZ?) into M right after (a,/?) and insert f(a, ab) into F. 
ii: end if 

12: Remove (a,b) from the list M. 
13: end while 



The only subtlety of this algorithm is updating the list M. We want to remove (a,b), of 
which we know that it is a minimal element of Xj. The set of elements strictly greater than 
(a,b) is given by 

{{a,b)Y^{(a,b)} = {(a,ab),{aa,b)Y. 

Hence, removing (a,b) creates at most two new minima, (a,ab) and (oa,b). It is easy 
to check whether they are actually minimal for Xj \ {(a,b)}: since our list M of minima 
is ordered, it suffices to compare (a,(jb) to the predecessor (#*,&*), and (oa,b) to the 
successor (a*,b*). 

To illustrate the different possibilities, we consider Figure [2] again. The following table 
indicates, for each possible minimum (a,b), how the list M has to be modified in order to 
obtain a new ordered list of minima satisfying M # — X, \ {(a,b)}: 
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yes 


(9,2) 


(12,2) 


no 



Lemma 12. Algorithm$\is correct: if A and B are finite ordered sets and f : Ax B — > Z is a 
bimonotone map, then Algorithm^produces a stream enumerating all pairs {a,b) G A x B 
such that the values f{a 1 b) appear in increasing order. 

Proof. At the beginning of the i-th iteration of the algorithm we denote M by M,-, and 
F by Fi, and the set of remaining parameters by Xj := Mf. The initialization states that 

Ml = ( (flmin, Vin) ) and F\ = ( /(fl min , Vin) ). SO X { = A X B. 

By induction we can assume that the set X,- is of size n — i+1 and saturated, with 
Mi = MinpQ) and F; = /(M,-). Furthermore, we can assume that the list representing M, is 
ordered in the sense that ( (a\ ,b\), (02,^2), ■ ■ ■ , (a m ,b m ) ) satisfies ai < a-i < ■■■ < a m and 
b\ > b% > ■■■ > b m . It is straightforward to verify that the 2-th iteration of our algorithm 
ensures the following assertions: 

• The output Xj satisfies /(x,) = minf; = min/(M,) = min/pQ). 

• The set Xj + i — X, \ {x,} is saturated and of size n — i. 

• We have M i+ \ = Min(Zj +1 ) and F /+1 = f(M i+1 ). 

• The new list representing the set M, + i is again ordered. 

The algorithm stops after « iterations when it reaches X n+ i — 0, hence M n+ i = and 
Fn+i = 0- We conclude that the output sequence Xi,X2, ■ ■ ■ ,x n is an enumeration of Ax B 
satisfying/^) < f(x 2 ) < < f(x n ). □ 

Since the sorted enumeration algorithm outputs one element x, at each iteration, the 
loop is repeated exactly n — \A \ ■ \B\ times. At the 2-th iteration, the algorithm occupies 
memory of size m; = |M, | to store the list M; and the priority queue F'. Let m = maxm; be 
the maximum during the entire execution. 

Lemma 13. In order to enumerate a set X = A x B of size n, Algorithm^\only builds up a 
priority queue of length m < min{|A|, \B\}, which entails in particular m < y/n. 

Proof. During the enumeration algorithm the list representing M is always strictly increas- 
ing in a and strictly decreasing in b. In particular, the projections M — > A and M — > B are 
both injective. The required memory m is thus bounded by min{ |A| , |5| }. □ 

Proposition 14. Let f: NxN->Zka bimonotone map that behaves polynomially. Ap- 
plied to subsets A = { 1 , . . . , 1} and B = { 1 , . . . , m} with m > 2, Algorithm [4] takes time 
O(n\og 2 n) and requires memory 0(^/nlogn). 

Proof. The loop is repeated n times. The most time consuming operation is updating 
the priority queue F; to F + i which requires time O (log n log m,), where m,- is the size 
of the queue F; an d its elements are typically of size ©(log 71). The total cost is time 
0(/i log n log ra) and memory 0(77tlogn). With m < y/n we obtain the stated bounds. □ 



10 



MICHAEL EISERMANN 



Example 15. The following extreme cases illustrate Algorithm @] and the possible be- 
haviour of the memory bound m. We consider /: A x B — > N where A = {l,...,k} and 
B = {1, ...,/} are two intervals of integers, with k > I > 2 say. 

The best case occurs for f(a,b) = la + b, where Algorithm |4] enumerates A x B in 
lexicographic order. During the z'-th iteration of the algorithm the set of minima M, contains 
only 1 or 2 elements, so that m = 2, independent of the sizes \A\ and \B\. 

The worst case occurs for f(a,b) = a + b. Having enumerated all elements x, with 
f(xi) < I, the list M contains exactly I elements, namely (1,1), (2, 1 — 1), ... , (/, 1). Thus 
the upper bound m = min{|A|, is actually attained. 

4.4. Enumerating infinite sets. As with semimonotone enumeration, bimonotone enu- 
meration can be generalized from finite to infinite sets. The interesting point is that now 
both sets A and B can be infinite, and the algorithm applies without change. 

Theorem 16. Suppose that A and B are ordered sets and that f : A x B — > Z is a proper 
bimonotone map. Then Algorithm\4\produces a stream enumerating all pairs (a,b) £AxB 
such that the values f{a,b) appear in increasing order. In order to enumerate the first n 
values, the algorithm builds up a priority queue of length at most \ 2n-\- 1. If f behaves 
polynomially, the algorithm takes time <9(«log 2 n) and requires memory 0(^/nlogn). 

Proof. For every z £ Im(/) the set { (a,b) £ A x B | f(a,b) < z} is finite, it is thus contained 
in some finite product [^min? a\] x [bmimbi]. Hence Algorithm Incorrectly enumerates all 
parameters (a,b) with f(a,b) < z as in the finite case. Since this is true for all z, the 
enumeration exhausts Ax B. 

Let us suppose that after n outputs the list M holds m pairs (ai,bi) ordered such that 
a\ < fl2 < ■ ■ ■ < fl;n and b\ > b% > ■■• > b m . This obviously implies n > lm(m— 1), whence 
m < V2n + 1. The time needed for n outputs is thus (9(nlog«log ^/n) = <9(nlog 2 n), while 
the required memory is O ( \fn log n ) . □ 

Example 17. Again the worst case occurs for the map / : NxN^Z with f(a,b) =a + b, 
where m ~ \/2n. The best case occurs for f(a,b) = max{a,b}, where m < 4. 

This example shows that for bimonotone enumeration the memory bound m £ 0(s/n) is 
best possible: there exist bimonotone functions / for which Algorithm [4] actually requires 
temporary memory m ~ y2n. Notwithstanding, the algorithm performs significantly better 
on certain subclasses of bimonotone functions: 

Proposition 18 (separate variables). Consider f : NxN^Z with f(a,b) = p(a) +q(b), 
where p and q are non-decreasing polynomial functions of degree a — deg p and j3 = 
degq, respectively. Assuming 1 < a < j3, the bimonotone enumeration algorithm requires 
memory m £ 0(n £ ) with exponent £ = ^nj. 

For example, sorted enumeration of f{a,b) = a 3 + b 1 requires memory m £ ©(n 3 / 10 ), 
while the a priori bound of Theorem[T6lonlv tells us m £ 0(n 1 ' 2 ). 

Proof. Let /: N x N -> N be defined by f(a,b) =a a +b^ with 1 < a < J3. The general 
case f(a,b) — p(a) +q(b) works essentially the same, but notation is more cumbersome. 

Suppose that the n-th output x„ has attained the level f(x„) — z, and the list M holds 
m parameters (a\,b\), . . . , (a m ,b m ). Then we have a\ = and f(a\,b\) = b\ > z. On the 
other hand {a\,b\ — 1) has already been output, which means (b\ — \)^ < z. We conclude 
that < b\ < ^/z+ 1, whence b\ ~ %fz. Analogously ^/z < a m < yfz+ 1, whence 
a m ~ \/z- This situation is depicted in Figure [3] In the sequel we set a : = a m and b:=b\. 
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FIGURE 3. Estimating the size of the set Mm^{a a +b& > z} 



The upper bound n < (a + 1) (b + 1) is clear. Since / is convex, we also have the 
lower bound n > jab. To see this, apply Pick's theorem to count integer points in the 
triangle A = [(0,0), (a,0), (0,b)]. Both inequalities together imply that n G ©( z ( a +P)/ a 0), 
or equivalently, z € &(n a ^^ a+ ^). 

The upper bound m < b + 1 is clear, and it remains to establish a lower bound. We 
will assume a < j3. (The symmetric case a = j3 is easier and will be examined more 
closely in Example |271 below.) There exists a unique point (x,y) £ on the contour 
whose normal vector points in the direction (1,1): this is the solution of x a + yP =z and 
ax a ~ l = jiyP~ l - It is easy to see that m>b— \y\ and m ~ (a — x) + (b — y) as indicated in 
Figure|3] We havex = cy^- 1 )/( a - 1 ) withe = a -{ffja, whence cV^- 1 '^"- 1 ^/ =z. 
Since a(j3 - l)/(a - 1) > j3 we deduce that y € o{ %/z). The bounds b-\y \ <m<b+l 
thus entail m — %/z, whence m G 0(n a /( a+ ' 3 )). □ 

We remark that in the above examples semimonotone enumeration achieves the same 
asymptotic bounds. This warrants a more detailed analysis, which we endeavour next. 

5. Asymptotic complexity 

We are now ready to address the crucial question: is bimonotone enumeration (Algo- 
rithm|4]i better than semimonotone enumeration (Algorithm|3]l? We shall compare the size 
m of the priority queue built up during the algorithm. The test class consists of all proper 
bimonotone functions / : N x N — > Q, which is where both algorithms apply. First of all, 
the following observation is worth emphasizing: 

Remark 19. Bimonotone enumeration is at least as good as semimonotone enumeration. 
More explicitly, both algorithms have to trace the contour of the finite set 

{f^z}:={(a,b)eAxB\f(a,b)^z} 

and construct the set of minima of the complement, Min{/ > z}- To this end semimonotone 
enumeration uses the partial order (a,b) < (a' ,b') defined by a ^ a 1 and b = b' . (Here we 
are ordering with respect to a for fixed b; since / is bimonotone, we could also order with 
respect to b for fixed a, whichever is advantageous.) Bimonotone enumeration uses the 
partial order (a,b) << (a' ,b') defined by a ^ a' and b ^ b'. This entails the inclusion 

Min<K{/ > z} C Min<{/ > z}. 
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This means that the priority queue for bimonotone enumeration is a subset of the queue for 
semimonotone enumeration, and consequently the required memory is less or equal. 

At this point we should clarify a possible ambiguity. Both Algorithms [3] and [4] have 
to choose one minimal element of the priority queue. In order to disambiguate multiple 
minima, we choose the one with minimal B-coordinate. This ensures that it belongs to both 
Min«{/ > z} and Min<{/ > z}, and the inclusion propagates inductively. 

Whether the bimonotone algorithm can achieve a significant improvement depends on 
the function /. Let us begin with a trivial example where no savings are possible: 

Example 20 (linear contour). Consider /: NxN->N defined by f(a,b) = (a + b) r with 
7 > 1 . In this case Min« {/ > z} = Min< {/ > z} is given by the line a + b = 1 + [-^/z\ . 

In general, however, the inclusion Min«{/ > z} C Min<{/ > z} is strict. Generally 
speaking, bimonotone enumeration adapts better to the contour and achieves savings when- 
ever the contour deviates from being a straight line. We now quantify this observation. 

5.1. Polynomial functions. Consider /: NxN->Q defined by a polynomial f{a,b) = 
Y.i,j c ij a 'bj with rational coefficients c, ; > 0. This condition ensures that / is bimonotone. 

Let fi(a) — f(a,0) and /2(b) — f{0,b) be the induced polynomial functions on the 
axes, and set a := deg/i and /3 := deg/2. We assume that a,/3 > 1, which ensures that / 
is proper. Without loss of generality we can also assume that a < /3 . 

Let 7 := deg/ = max{i + j | c,/ ^ 0} be the total degree of /. We have a < j5 < 7. 

We denote by n := ^ z} the number of values of / up to some level z, and by 
m := jjMin{/ > z} the length of the priority queue at level z- 

Proposition 21. Semimonotone enumeration of the set {f ^ z} requires memory m £ 
0( ^/z) whereas bimonotone enumeration requires memory m € 0{^/z). 

Whenever j3 < 7, bimonotone enumeration is thus significantly better than semimono- 
tone enumeration. As an example, for f(a,b) ~a 4 + a^b 4 + b 5 semimonotone enumeration 
requires memory m 6 &(^/z), while bimonotone enumeration requires only m E 0(J/z). 

Proof. Assuming j8 > a, it is advantageous to sort by a in the semimonotone enumeration 
algorithm. In this case we see that m = 2 + b where b satisfies /2(b) < Z < fiip + l). We 
have fi (b) ~ cb@ with some leading coefficient c > 0, whence m ~ {/z/c. 





b 




x 



X 



a 



Figure 4. Estimating the size of the set Min«{/ > z} 
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Evaluating / on the diagonal, we find f(x,x) = do + d\X H h dyX^ with non-negative 

coefficients dk>0 and dy > 0. For x = 1 + [{/z/dy\ we obtain f(x,x) > z. This is illus- 
trated in Figure |U where the dotted line corresponds to / = z, and black dots represent the 
elements of Min«{/ > z}. We conclude that m < 1 + 2x, whence m 6 O(-^z). □ 

5.2. Asymptotic bounds. In order to express the required memory m in terms of the 
number n of enumerated values, we wish to relate n and z. For z — > 00 we can replace 
counting points (a,/?) G N 2 satisfying f(a,b) < z by the Lebesgue measure of the set 

{(x,y)eR 2 + \f(x,y)<z}. 

Proposition 22. Let f: — > R + /?e a polynomial function given by 

f(x,y) = c ij x 'y J with Cjj > for all indices (i,j) G K. 

(i,j)eK 

With f we associate the convex polygon D = {(a,v)£l 2 | iu + yv < 1 /or aZZ (i, j) G A'}. 
Suppose that f is proper in the sense that for all z G M+ the set 

{f<z} = {(x,y)eR 2 + \f(x,y)<z} 

is compact. Then its Lebesgue measure satisfies A({/ < z}) G &(z S log(z) d ) where 

8 := max{(( + v | (m, v) G £>} 

anrf J is the dimension of the set where this maximum is attained: either d = Ofor a vertex, 
or d = 1 for a segment. (See Figure\5\belowfor examples.) 

The proof is a nice application of the so-called "tropical" approach. The idea is to 
identify R+ = {x G K | x > 0} and R = R U { -°°} via the natural logarithm log : R+ -»■ R, 
and to formally replace the semiring (R+, +, •) by the semiring (R.max, +). Of course, we 
have \og(x-y) — logx + logy but for log(x + y) we only obtain an inequality, 

max(logx.logy) < log(x + y) < log2 + max(logx,logy). 

This means that log : R + — > R is a quasi-isomorphism, i.e. its failure to be an isomorphism 
is bounded by some constant. For asymptotic arguments this is usually sufficient. 

Proof of Proposition^^ As a logarithmic analogue of / we define 

/: R 2 - R, f(x,f) := max (ix + jy) . 

We can choose a constant c G R+ such that cy > e~ c and (tl^T) • cy < e +c for all (z", j) G A". 
A small calculation then shows that 

|log/(x,y)-/(logjc,logy)| <c. 

The measure of the set {/ < z} equals the integral over the associated indicator function 
[/ < z] and we can apply the change of variables x — logx, y — logy, z — logz: 

F (z) ■= / , if(x,y) < z] dxdy = / [log/(i,y) < z]exp(x + y)dxdy. 

Jr\ Jr 2 

It is easier to calculate this integral with / instead of /, so let us do this first. Since / is 
homogeneous, we perform another change of variables x = uz and y = vz to obtain: 

F(z):= / \f(x,y) < zl exp(i + y)dxdy = z 2 / [/(a, v) < ll exp(wz + vz)dudv. 
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We are now integrating over the convex polygon D := {(u,v) £ R 2 | f{u,v) < 1}. The 
asymptotic behaviour of logF(z) is easy to understand: for z — ► °° we obtain 



logz 



log 



?2/z 



exp(t< + v) dudv 



For z — > °° the first factor z 2// ~ — ► 1 plays no role. The remaining factor is the z-norm 
|| exp(w +v)||| and tends to the sup-norm || exp(u + v) ||oo = exp(5) for z — ► 
This shows that \ogF{z) ~ 5 logz, but does not yet suffice to imply F (z) 
We thus have a closer look at the quotient 

Hz) 



<z S forz 



log(z) 2 / z {u+v - S) dudv. 
Jd 



We change variables u 



S-t 
2 



s and v : 



5-f 

2 



s so that k - 



8 — t and dudv = dsdt: 



> +v - 5 Uudv= / l{t)z~*dt 
id Jo 

where t(t) is the length of the segment {(«, v) G D | u + v = 8 — t }. Since D is a polygon, 
there exist ao,fli > and T > such that ^(f ) = aq + «if for all t £ [0, T]. We thus find 



/ £(t)z 'dt~a log(z) 
Jo 



-l 



-ailog(z) 



-2 



forz 



Notice that «o = if and only if the maximum u + v = 8 is attained in a single vertex. 
We thus obtain F(z) £ &(z s log(z) rf ) where 5 is the maximum of u + v on D and d is 
the dimension of the maximising set. Since F(e~ c z) < F(z) < F{e +C z), we conclude that 
F{z)£&{z s log(z) d ). " □ 

Remark 23. It is clear that the proposition and its proof generalize to proper polynomial 
functions /: M'| — > R + with non-negative coefficients, in any number n of variables. We 
have concentrated on n = 2, which is the case of interest to us here. 



Example 24. For f(x,y) 



-y 5 the set {/< 1} is depicted in Figure|5]on the left. Here 



we obtain 8 = and d = because the maximum is attained in a single vertex. 

The figure in the middle shows {/ < 1} for f(a,b) = a 4 +a 3 fr 3 +b 5 . Here 8 = I and 
d = 1, because the maximum is attained on a segment, so that n £ @(z 1,/3 logz). 



(1/4,1/5) 



(2/15,1/5) 




1/4,1/12) 




(1/4,1/16) 



Figure 5 . Maximizing u + v under the constraint f(u,v) < 1 

The figure on the right shows {/ < 1 } for f(a,b) = a 4 + a^b 4 + b 5 . Here we find 8 = , 
which means that z £ @(« 16 / 5 ). According to Propositionl2"Tl semimonotone enumeration 
requires memory m £ ©(z 1 / 5 ), whence m £ @(« 16 / 25 ). 

Notice in particular that > \- This illustrates that, unlike bimonotone enumeration, 
semimonotone enumeration cannot guarantee the memory bound m £ Oi^fn). 
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Corollary 25. The semimonotone enumeration algorithm guarantees the memory bound 
m<n + \, and as a uniform bound the exponent 1 is best possible: enumerating the values 
of f{a,b) = a a +a a bP +b^ with a < j3 requires memory m £ &(n a ^). □ 

Proof. According to Proposition [22] the number of enumerated values up to level z is 
n £ &(z S ) with 8 = l/a, and thus z £ ©(«"). According to Proposition |2T1 the required 
memory is m £ &(z 1 ^). We conclude that m £ ®(n a 'P). □ 

5.3. Constant factors. Propositionl2Tlexhibits many polynomial functions where bimono- 
tone enumeration is clearly worth the effort. Depending on the envisaged application and 
the given function /, a finer analysis and a more modest conclusion may be necessary: 

Example 26. Consider polynomials of the form f(a,b) = p(a) +q(b), for which semi- 
monotone enumeration was initially devised H][Q. We obtain n £ &(z S ) with 8 = ^ + 4, 
as already remarked in the proof of Proposition [18] Assuming a < j3, bimonotone and 
semimonotone enumeration both require memory m £ 0(n £ ) with £ = ^^j- 

Even if memory requirements are of the same order of magnitude, we can usually expect 
to gain a constant factor with the bimonotone algorithm: 

Example 27. Reconsider /: NxN->Z defined by f{a,b) =a r + b r with y > 1. In this 
case semimonotone enumeration requires memory m ~ \fz, whereas bimonotone enumer- 
ation requires memory m ~ cy tfz with a factor cy = 2(1 — %J 1 /2) < 1 . 

Proof. For a = b = L^/zJ we have f\ (a) = fi{b) < z and fi(a + 1) = fi(b +1) > z, and 
semimonotone enumeration requires memory (JMin<{/ > z} = 2 + b ~ ^fz. 

Choosing x = y = [-^/z/2\ we find f(x,y) < z and /(x+l,y+l) > z. As indicated in 
Figure[3] we have m ~ (a — x) + (b — y) ~ Cy ■ ^fz. □ 

Though less impressive, for practical applications even a constant factor may be a wel- 
come improvement: reducing memory consumption means that we can scale to consid- 
erably larger problems before running out of RAM. In our example we have C2 ~ 0.59, 
c 3 ps 0.41, c 4 pa 0.32, c 5 ps 0.26, and c y -> for y ^ °°. 

6. Parallelization 

Let us reconsider the application of sorted enumeration to a diophantine equation f(a,b) = 
g(c,d), where f,g: NxN->Z are proper bimonotone functions. Suppose we are look- 
ing for solutions x — (a,b), y — (c,d) with values in some large interval z m ; n < f(x) = 
g(y) < Zmax- This problem can be split into s independent subproblems, namely search- 
ing solutions with Zk-i < f(x) = g{y) < Zk, where z min = zo < z\ < zi < ■ ■ ■ < z s = z ma x 
is a subdivision of our search interval. This allows us to distribute the search on several 
computers in parallel. 

6.1. The initialization algorithm. To put the parallelization idea into practice, Algorithm 
[5] stated below, initializes the enumeration stream to begin at level z. Graphically speaking, 
it traces the contour of X = {/ > z} in order to determine the set of minimal elements 
M = MinX. From M we can then immediately build up the priority queue F = f(M). 

As usual we require that /: A x B — > Z be a proper bimonotone map. For simplicity we 
first assume that both A and B are infinite. (We will treat the general case in the next para- 
graph.) As before the successor function is denoted by a 1— > aa and b ab, respectively. 
We also use the predecessor function, denoted by £» 1 — %b. 
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Algorithm 5 Constructing the set of minima on X = A x B 
Requires: a proper bimonotone function / : A x B — * Z 
Input: a level z £ Z 

Output: the list of minima Min{x £ X \ f(x) ^ z} = ( {a\,b\), (02,^2), ■ ■ ■ , (a m ,b m ) ) 
Ensures: a\ < «2 < ■ ■ ■ < a m and b\ > bi > ■ ■ ■ > b m 

l: Initialize M <— and a <— a m i n , b <— b m ; n 

2: while f(a,b)<z do b^ab end while 

3: Insert (a,&) at the end of the list M; continue with a <— era, b ^ Kb 
4: while > fr min do 

5: while f(a,b) <z do a<— aa end while 

6: while > fomin and f(a, Kb) >z do b Kb end while 

7: Insert (a,b) at the end of the list M; continue with a <— era, b ^ Kb 

8: end while 

9: return M 



The reader is invited to apply Algorithm|5]to the example given in Figure|2] in order to 
see how it traces the contour of X = {/ > z}- By the way, the method applies to any set 
X C Ax B that is saturated and has finite complement. We shall give a detailed proof in 
the more general situation of Algorithm |7]below. 

Remark 28. The loop in line 2 determines b «— min{Z? £ B | f(a m i n ,b) ^ z}. This could, 
of course, be improved by replacing the linear search with a binary search, provided that 
b can easily be incremented and decremented by integer values. The same holds for the 
loops in lines 5 and 6. This optimization is straightforward to implement whenever the 
application requires it. 

Remark 29. Let M = ( (ai,b\), (02,^2), ■ ■ ■ > (am,b m ) ) be the list of minima. Then build- 
ing a priority queue from M requires time O(mlogn). Moreover, let k = #[a m i n ,a m ] and 
l=#\b min:^i]? with k > /, say. Then n > k > / > ?w. Constructing the list itself requires 
time (9(fclogn) using linear search, and 0(ralog 2 n) using binary search. We cannot expect 
to do much better, because constructing a list of length m requires at least m iterations. 

6.2. Applications. Having initialized M and F, we can apply the bimonotone enumeration 
algorithm to produce a sorted enumeration x\,X2, ... of the set {/ > Zk-i}- Applying the 
same method to g, we can produce a sorted enumeration yi,y2, ■ ■ ■ of {g > Zk-i}- We can 
thus search for solutions f(x) = g(y) starting at level Zk-\ and ending at level Zk- 

Expected speed-up. Concerning time requirements, initialization entails a reasonably 
small overhead, so we can expect an amortized speed-up by a factor s. For each k = l,...,s, 
computer number k manages its own priority queues of length 0(y/n). in order to produce 
enumeration streams for / and g, with values ranging from Zk- 1 to Zk- As before, advancing 
from position n to position n + 1 takes time <9(log 2 n). 

Robustness. The initialization procedure is already very useful on a single computer, 
since it can make implementations much more robust: it is possible to continue searching, 
without much loss, after a shut-down or a power failure. This is particularly important 
when carrying out a long-term search. 

7. Enumerating bimonotone domains 

Suppose we want to enumerate the values of a symmetric bimonotone function / : N x 
N — > Z, that is, f(a,b) = f(b,a) for all a,b £ N. It is often desirable to enumerate only 
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pairs (a,b) with a^b. In other words, we wish to restrict f to the domain X = { (a,b) £ 
N x N | a ^ b } and enumerate only the values of / : X — ► Z. Figure [6] shows a possible 
configuration during bimonotone enumeration. 
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FIGURE 6. Enumerating the domain X = {(a,Z?)eNxN|a^fr} 

It is straightforward to adapt bimonotone enumeration (Algorithm 2]i and initialization 
(Algorithm [5]) to such a domain X. Still other restrictions are possible, for example X = 
{ (a, b) E N x N \ a < b and b <2a} or more complicated cases such as X = { (a,b) £ 
N x N | a < b and b 2 < 1 + 10a }. This raises the question as to what are "reasonable" 
domains X C A x B to which Algorithms[4]and[5]can be efficiently applied. 

7.1. Bimonotone domains. As usual we assume that A and B are isotonic to finite inter- 
vals or to the natural numbers. Figure [7] shows a domain X cA x B which will turn out 
to be well suited to bimonotone enumeration. Graphically speaking, it is bounded by the 
graphs of two non-decreasing functions a : A — > B and j3 : B — > A. We will show that this 
condition suffices to adapt our algorithms to work on the domain X rather than the entire 
product A x B. 



| X X X X X X 

fx X X X X X X 
XXXXXXXX 

xxxxxxx 

® X X X X X x 

• X X X X X 

• ® X X X 



H 1 h 



H h 



H 1 h 



10 11 12 13 14 15 



Figure 7 . Enumerating a bimonotone domain X cAx B 

We say that X C A x B is bounded by functions a : A — > B and /3 : B — > A if 
X = { (a,b) e A x B | a > ]8 (b) and b > a {a) }. 
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For example, the domain X of Figure[7]is bounded by a(l) = ■ ■ ■ = a(12) = 1, a(13) = 4, 
a(14) = 5, a(15) = 7, and j8(l) = 1, 0(2) = 0(3) =3 J3(9) = 10. 

Definition 30. We say that a domain X C A x B is bimonotone if it is bounded by two 
functions a.A^B and fi.B^A such that: 

(1) The functions a and are non-decreasing, that is, 

a < a' implies a(a) < a(a'), and b < implies 0(7?) < J3(fo'). 

(2) We have 0(a(a)) ^ a for all a S A, with equality only for a = a m \ a , 
and a(P(b)) ^b for all b S B, with equality only for = £> m j n . 

Condition (2) ensures that (a,a(a)) £ X for each a G A, and (f}(b),b) £ X for each 
€ 5. In particular, a and are determined by X via 

a(a) = min{Z? e B \ (a,b) £X], 
P(b)=min{a£A\{a,b)£X}. 

Moreover, (a m i n ,£> m i n ) is the smallest element of X. If both A and B are finite, then 
(flmax,^max) is the greatest element of X. 

The definition of X via bounding functions is easy to formulate and well suited to im- 
plementation. It can also be reformulated in more geometric terms: 

Proposition 31. A domain X C Ax B is bimonotone if and only if it satisfies pi"j X — A and 
pr 2 X = B and the following two properties: 

(1') If (a\,b\) and (a^bi) in X satisfy a\ ^a2andb2 ^b\, then X contains the entire 

rectangle { (a,b) £ A x B \ a\ ^ a «2 andbi ^ b ^ b\}. 
(2') 7f ant/ {0,2^2) in X satisfy a\ «2 ant/ £>i $J £>2> we can go from 

{a\,b2) to (ai^bj) within X by repeatedly incrementing a and b. □ 

The proof is not difficult and will be omitted. 

7.2. Bimonotone enumeration. We are now in position to generalize our enumeration 
algorithm to a bimonotone domain. As before, Algorithm|6]processes a bidirectional list 
M and a priority queue F. 

Proposition 32. Algorithm]6\is correct. 

Proof. The proof is essentially the same as for Algorithm |4] There are, however, some 
modifications when updating the list M and the priority queue F: 

• If the current minimum (a,b) is somewhere in the middle of the list M, then the 
previous arguments apply without change, because we still have 

{(a,b)} # \ {(a,b)} = { (a, ob), (oa,b) } # . 

• If (a,b) is at the end of the list, then possibly (oa,b) ^ X: in this case we have 
{(a,b)} # \ {(a,b)} — {(a,ob)} # , so we discard (oa,b). 

• If (a,b) is at the beginning of the list, then possibly (a,ob) £ X: in this case we 
have {(a,b)} # \ {(a,b)} — {(oa,b)} # , so we discard (a, ob). 

• If ever M — ( (a,b) ) and neither (a,ob) nor (oa,b) is in X, then (a,b) is the 
greatest element of X and the algorithm terminates correctly. 

Since / is proper, every element (a,b) £ X will eventually be enumerated. □ 
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Algorithm 6 Sorted enumeration of a bimonotone domain 

Requires: a bimonotone domain X cAxB and a proper bimonotone function / : X — ► Z 
Output: an enumeration of X, monotone with respect to / 

1: Initialize M <— ( {a min7 b min ) ) and F <— ( /(a m i n ,&min) )■ 

2: while F is non-empty do 

3: Remove a minimal element f(a,b) from F and output (a, ft). 

4: if (a, fo) is the last element of the list M then 

5: if (oa,b)&X then insert [aa,b) into M and f(oa,b) into F end if 

6: else 

7: Let (a* ,b*) be the successor of (a,£>) in the list M. 

8: if oa<a* then insert (oa, b) into M and /(era, fe) into F end if 

9: end if 

10: if {a, b) is the first element of the list M then 

11: if (a, afo) £ X then insert (a, ob) into M and /(a, ob) into F end if 

12: else 

13: Let {a*,b*) be the predecessor of (a,b) in the list M. 

14: if ob < b* then insert (a, a/?) into M and f(a, ob) into F end if 

15: end if 

16: Remove (a,b) from the list M. 

17: end while 



7.3. Initialization on a bimonotone domain. As for the unrestricted case X =Ax B,we 
want to formulate an initialization algorithm for a proper bimonotone function f:X—> 
Z defined on some bimonotone domain X C A x B. The idea is essentially the same: 
Algorithm [7] traces the contour of X(z) = {x € X | f(x) ^ z} to construct the list M = 
MinX(z) of its minimal elements. 



Algorithm 7 Constructing the set of minima on a bimonotone domain 

Requires: a bimonotone domain X cAxB and a proper bimonotone function / : X — > Z 

Input: a level z € Z 

Output: the list of minima Min{x € X | /(*) ^ z} = ( (02,^2), ■ ■ ■ > (a>n,b m ) ) 

Ensures: «i < «2 < ■ ■ ■ < &m an d b\> bi> ■ ■ ■ > b m 

1: Initialize M <— and a <— a m i n , <— b m { n 

2: while (fl,i)eA r and/(fl,l))<z do 

3: while f(a,b) < z and (a, ob) EX do b <— ob end while 

4: if f(a,b) < z then oa end if 

5: end while 

6: while (a,!))6X do 

7: while (a, Kb) e X and f(a, nb) ^ z do b^%b end while 

8: Insert (a,b) at the end of the list M; continue with a «— aa, b ^ nb 

9: while (a,fc)GXand/(a,i?)<z do a^oa end while 

10: end while 

li: return M 



Remark 33. The loops in lines 3, 7, and 9 implement linear searches. This can be im- 
proved by a binary search whenever the application requires such optimization. 

Proposition 34. Algorithm\7\is correct. 
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Proof. The first loop (lines 1-5) finds an element (a,b) £ X(z) with minimal a. Beginning 
with a «— fl m i n and b *— b m \ B , we repeatedly increment b in order to arrive at f(a,b) z. 
If this is not possible within X, then the candidate a is eliminated, and we continue with 
a «— <Ja. If we never run out of the domain X, then we finally end up with f(a,b) ^ z, 
because a or b increase and / is proper. 

The only obstacle occurs when f(a,b) < z but neither (a, ab) nor (aa,b) are in X. But 
in this case we have reached the greatest element of X, hence f(x) < z for all x £ X. Thus 
X(z)=% and we correctly return the empty list M= 0. In any case, the first loop terminates 
with either (a,b) £ X or f(a,b) ^ z, as desired. 

When arriving at line 7 we know that (a, b) £ X(z) \M # , and a is minimal with this 
property. The loop in line 7 minimizes b, so we know that (a,b) is a minimal element of 
X(z). We thus add (a,b) to our list M and continue with a «— era and Z? <— 7r£>. We then 
repeatedly increment a in order to arrive at f(a,b) ^ z- If this is not possible in X, then 
X(z) — M* by the rectangle condition (T), so we have found all minimal elements of X(z). 
Otherwise, we obtain (a,b) £ X(z) \M # , and a is minimal with this property. We can thus 
reiterate by looping back to line 7. 

During each iteration, a is strictly increasing while b is strictly decreasing. We conclude 
that the second loop terminates and produces the list M of minima, as desired, ordered in 
the sense that a\ < 02 < ■ ■ ■ < a m and b\ > bi > ■ ■ ■ > b m . □ 



Algorithms [6] and [7] for bimonotone enumeration have been implemented as a class 
template in C++. This seems to be a good compromise between general applicability, ease 
of use, and high performance. The source files are available on the author's homepage: 
http : / / www-f ourier .ujf-grenoble. f r / ~eiserra/ software 

As an illustration of sorted enumeration, let us mention searching multiple values of a 
polynomial function /: N x N — * Z, f(a,b) = Y^i^cijdV with non-negative coefficients 
Cij £ N. The cited implementation has been successfully tested to reproduce some known 
results taken from Richard Guy's Unsolved problems in number theory [4]. 

8.1. The quest for the sixth taxicab number. As an illustrative example we briefly 
sketch the taxicab problem. The kth taxicab number, denoted by taxicab(fc), is the least 
positive integer that can be expressed as a sum of two positive cubes in k distinct ways, up 
to order of summands. That is, it is the smallest £-fold value of f(a,b) = a 3 +b 3 defined 
onX = {(fl,i)eNxN 1 <a<b}. 

G.H.Hardy and E.M.Wright [5 Thm.412] proved that, for every k > 1, there exist 
such £-fold values. This guarantees the existence of a least £-fold value, that is, the &th 
taxicab number. Unfortunately the construction given in the proof is of no help in finding 
the least fc-fold value. Apart from (variants of) exhaustive search, no such method is known 
today. The first taxicab number is trivially 



8. Applications to diophantine enumeration 



taxicab(l) 



2 = 1 3 + 1 3 . 



The next taxicab numbers are: 



taxicab(2) = 1729 



1 3 + 12 3 =9 3 + 10 3 
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(re)discovered by Ramanujan according to Hardy's famous anecdote, but previously pub- 
lished by Bernard Frenicle de Bessy in 1657, 

taxicab(3) = 87539319 

= 167 3 +436 3 = 228 3 +423 3 = 255 3 +414 3 , 

discovered by John Leech fH in 1957, 

taxicab(4) = 6963472309248 = 2421 3 + 19083 3 = 5436 3 + 18948 3 

= 10200 3 + 18072 3 = 13322 3 + 16630 3 , 

discovered by E. Rosenstiel, J. A. Dardis, and C.R. Rosenstiel f9) in 1991, 

taxicab(5) = 48988659276962496= 38787 3 + 365757 3 = 107839 3 + 362753 3 

= 205292 3 + 342952 3 = 221424 3 + 336588 3 = 231518 3 + 331954 3 , 

discovered independently by D . W. Wilson [ 1 1 in 1 997 and shortly afterwards by D . J . Bernstein 
Q] in 1998. Finally, the smallest known 6-fold value is 

T = 24153319581254312065344 

= 28906206 3 + 582162 3 = 28894803 3 + 3064173 3 = 28657487 s + 8519281 3 

= 27093208 3 + 16218068 3 = 26590452 3 + 17492496 3 = 26224366 3 + 18289922 3 , 

found by R.L. Rathbun in 2002. Is this actually the sixth taxicab number, or is there a 
smaller solution? 

8.2. Feasibility of an exhaustive search. In order to verify that T is indeed the smallest 
6-fold value, there are exactly n = 369039037733 393 < 4 ■ 10 14 pairs (a,b) G N x N to 
be checked with a 3 + £> 3 < T and a<b. Such counting results can easily be obtained from 
Algorithm|7]tracing the contour of the set X — {/ < z}: as a by-product, the initialization 
can be used to determine the sizes n = < z} and m — jjMin{/ > z}. 

Memory requirements are, fortunately, no problem. In the worst case we would have 
to check all n parameters, which would finally build up a priority queue of size m = 
5963 352 < 6 ■ 10 6 . Notice that each entry requires 32 bytes: 12 bytes for the value f(a,b), 
4 bytes for a and 4 bytes for b, plus 4 bytes for each of the three pointers. In the worst case 
the priority queue thus requires 180 megabytes of memory, which fits nicely in a PC with 
256 megabytes RAM. Such memory requirements seem acceptable; on today's PCs such a 
task can reasonably be run in the background. 

Time requirements, however, are on the edge of being feasible. Updating a priority 
queue of 2 • 10 6 entries, say, takes about 4000 CPU cycles. On a PC running at 2GHz, 
we can expect to process about 500000 steps per second, that is around 4 • 10 10 steps per 
day. This is not too far away from 4 • 10 14 , but on a single computer the search would 
still require about 10000 days, roughly 25 years. On 25 computers, however, we would be 
done within a year, possibly earlier. 

Partial results. Up to June 2005, I have run the search on a few available PC's at the 
Institut Fourier, but the use of parallelization has still been rather limited (to a dozen PC's). 
As a result I obtained the lower bound taxicab (6) > 5 • 10 20 by sorted enumeration of the 
2.8 • 10 13 smallest values of f(a,b) = a 3 +b 3 . (At a speed of 500000 values per second 
this takes about 650 days on a single computer.) This leaves us with the inequality 

5 ■ 10 20 < taxicab(6) < T w 2.42 ■ 10 22 

It will now be a matter of sufficient hardware and patience to find the exact answer. 
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